!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
!!
!! LDGH is free software: you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! LDGH is distributed in the hope that it will be useful, but WITHOUT
!! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
!! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
!! License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with LDGH. If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>


!> \brief
!! This is the classical cavity test case on the box \f$[0\,,1]^d\f$.
!!
!! \n
!!
!! The initial condition is identically zero; we enforce Dirichlet
!! boundary on all the faces of the box and it is assumed that the
!! boundary marker are set as follows:
!! <ul>
!!   <li> 1 for the upper lid, where we have
!!   \f$\underline{u}_D=U\underline{e}_1\f$, for some
!!   \f$U\in\mathbb{R}\f$ and with \f$\underline{e}_1\f$ denoting the
!!   unit vector along the first coordinate axis;
!!   <li> any marker \f$\geq2\f$ for the remaining boundary regions,
!!   where we set \f$\underline{u}_D\equiv0\f$.
!! </ul>
!!
!! \warning Since there are only Dirichlet boundary conditions, it is
!! advisable to use a FE method which is robust with respect to the
!! pressure computation and the incompressibility constraint.
!!
!! There is no known analytic solution for this test case; some
!! reference solutions can be found in <a
!! href="http://dx.doi.org/10.1016/0169-5983(89)90020-8">[Iwatsu,
!! Kawamura, Kuwahara, Hyun, 1989]</a>, <a
!! href="http://dx.doi.org/10.1016/S0045-7930(98)00002-4">[Botella,
!! Peyret, 1998]</a>, <a
!! href="http://dx.doi.org/10.1016/j.jcp.2004.12.024">[Albensoeder,
!! Kuhlmann, 2005]</a>, <a
!! href="http://dx.doi.org/10.1002/fld.953">[Erturk, Corke,
!! G&ouml;k&ccedil;&ouml;l, 2005]</a> and <a
!! href="http://dx.doi.org/10.1002/fld.1887">[Erturk, 2008]</a>.
!!
!! \section numres Some numerical results
!! We fix here \f$U=1\f$ and vary the Reynolds number by changing
!! \f$\nu\f$.
!!
!! \subsection case2D The 2D case
!! We consider here the two-dimansional case.
!!
!! \subsubsection Re1000 The case Re=1000
!!
!! A grid converged solution can be obtained with
!! \f$\mathbb{P}_3-\mathbb{P}_3\f$ elements with pressure
!! stabilization with \f$h\approx 0.003962\f$, so that the maximum
!! difference with the data reported in tables 9 and 10 of <a
!! href="http://dx.doi.org/10.1016/S0045-7930(98)00002-4">[Botella,
!! Peyret, 1998]</a> is 0.005, with the sign corresponding to a weaker
!! vortex. For \f$t\to\infty\f$, the solution converges to steady
!! state with
!! \f{displaymath}{
!!  f-f_\infty \approx Ce^{-\alpha t}
!! \f}
!! with \f$\alpha\approx0.06\f$. A solution with
!! \f{displaymath}{
!!  \frac{f(t_1)-f(t_2)}{t_1-t_2} \approx 10^{-6}
!! \f}
!! is obtained for \f$T^{fin}\approx200\f$.
!<----------------------------------------------------------------------
module mod_cavity_test

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_fu_manager, only: &
   mod_fu_manager_initialized, &
   new_file_unit

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_cavity_test_constructor, &
   mod_cavity_test_destructor,  &
   mod_cavity_test_initialized, &
   test_name, &
   test_description,&
   coeff_visc,&
   coeff_f,   &
   coeff_dir, &
   coeff_init

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members
 character(len=*), parameter ::    &
   test_name = "cavity"

! Module variables

 ! public members
 character(len=100), protected ::    &
   test_description(3)
 logical, protected ::               &
   mod_cavity_test_initialized = .false.
 ! private members
 real(wp) :: nu, ub ! viscosity and imposed velocity
 character(len=*), parameter :: &
   test_input_file_name = 'cavity_test.in'
 character(len=*), parameter :: &
   this_mod_name = 'mod_cavity_test'

 ! Input namelist
 namelist /input/ &
   nu, ub

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_cavity_test_constructor(d)
  integer, intent(in) :: d

  integer :: fu, ierr
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.) .or. &
       (mod_kinds_initialized.eqv..false.)    .or. &
       (mod_fu_manager_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_cavity_test_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   ! Read input file
   call new_file_unit(fu,ierr)
   open(fu,file=trim(test_input_file_name), &
      status='old',action='read',form='formatted',iostat=ierr)
    if(ierr.ne.0) call error(this_sub_name,this_mod_name, &
      'Problems opening the input file')
    read(fu,input)
   close(fu,iostat=ierr)

   write(test_description(1),'(a,i1)') &
     'Cavity test case, dimension d = ',d
   write(test_description(2),'(a,e9.3,a)') &
     '  viscosity nu = ',nu,';'
   write(test_description(3),'(a,e9.3,a)') &
     '  boundary velocity Ub = ',ub,'.'

   mod_cavity_test_initialized = .true.
 end subroutine mod_cavity_test_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_cavity_test_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_cavity_test_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_cavity_test_initialized = .false.
 end subroutine mod_cavity_test_destructor

!-----------------------------------------------------------------------
 
 pure function coeff_visc(x) result(nu_x)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: nu_x(size(x,2))

   nu_x = nu
 end function coeff_visc
 
!-----------------------------------------------------------------------
 
 pure function coeff_f(x) result(f)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: f(size(x,1),size(x,2))

   f = 0.0_wp
 end function coeff_f
 
!-----------------------------------------------------------------------

 pure function coeff_dir(x,breg) result(d)
  real(wp), intent(in) :: x(:,:)
  integer, intent(in) :: breg
  real(wp) :: d(size(x,1),size(x,2))

   d = 0.0_wp
   if(breg.eq.1) d(1,:) = ub
 
 end function coeff_dir
 
!-----------------------------------------------------------------------

 pure function coeff_init(x) result(u)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: u(size(x,1),size(x,2))

   u = 0.0_wp
 
 end function coeff_init
 
!-----------------------------------------------------------------------

end module mod_cavity_test

